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^BSTRACT: Systematic numerical simulations of model dense granular materials in monotonous, quasistatic 
i— deformation reveal the existence of two different regimes. In the first one, the macroscopic strains stem from the 
^deformation of contacts. The motion can be calculated by purely static means, without inertia, stress controlled 
CN)r strain rate controlled simulations yield identical smooth rheological curves for a same sample. In the second 

! regime, strains are essentially due to instabilities of the contact network, the approach to the limits of large 

,I^amples and of small strain rates is considerably slower and the material is more sensitive to perturbations. 

^Fhese results are discussed and related to experiments : measurements of elastic moduli with very small strain 

increments, and slow deformation (creep) under constant stress. 

73 : 



INTRODUCTION 



Respite its now widespread use ([Kishino 2001), dis- 
fete numerical simulation of granular materials, mo- 
fivated either by the investigation of small scale (close 
i — to the grain size) phenomena, or by the study of mi- 

croscopic origins of known macroscopic laws, still 

graces difficulties. Microscopic parameters, some of 
C^vhich are to be defined at the (even smaller) scale 
^pf the contact, are incompletely known. Macroscopic 
^-constitutive laws do not emerge easily out of noisy 
_£imulation curves, and the numerically observed dy- 
t^yiamic sequences of rearrangements might appear to 
Ot;ontradict the traditional macroscopic quasistatic as- 
sumption. Detailed and quantitative comparisons with 
^experiments can be used to adjust microscopic mod- 
k^ls, but a systematic exploration of the effect of the 
M/arious parameters throughout some admissible range 
^s also worthwhile. This is the purpose of the present 
study, which also addresses the fundamental issues 
of the macroscopic and quasistatic limits, in the case 
of the biaxial compression of dense, two-dimensional 
(2D) samples of disks. 

In section 2, we introduce the model and the nu- 
merical methods and define dimensionless parameters 
that are robust indicators of the relative importance 
of different phenomena. Rheological curves can be 
evaluated in the large sample limit (section 3), and 
their sensitivity to parameters assessed. We observe 
(section 4) two different mechanical regimes, accord- 
ing to whether the dominant microscopic origin of 



strain is material deformation in the contacts or re- 
arrangements of the contact network. Connections to 
some experimental observations are suggested in part 
5, while the conclusion section outlines further per- 
spectives. 

2 NUMERICAL MODEL AND PROCEDURES 

2.1 Grain-level mechanics 

Our computational procedure is one of the simplest 
types of 'molecular dynamics' or 'discrete element' 
method (Cunda ll and Strack 1979|) for solid grains. 
We consider 2D assemblies of disks, with diameters 
uniformly distributed between a/2 and a, and masses 
and moments of inertia evaluated accordingly (as for 
homogeneous solid cylinders of equal lengths), m 
will denote the mass of a disk of diameter a, and N 
the number of disks. 

These grains interact in their contacts with a linear 
elastic law and Coulomb friction. The normal con- 
tact force Fjsi is thus related to the normal deflec- 
tion (or apparent interpenetration) h of the contact as 
F N = K N hY(h), Y being the Heaviside step func- 
tion (equal to 1 for h > 0, to otherwise). The tangen- 
tial component F T of the contact force is proportional 
to the tangential elastic relative displacement, with a 
tangential stiffness coefficient Kt- The Coulomb con- 
dition | Ft | < //-Fjv requires an incremental evaluation 
of F T every time step, which leads to some amount of 
slip each time one of the equalities F T = ±fxF N is 
imposed. A normal viscous component opposing the 
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relative normal motion of any pair of grains in contact 
is also added to the elastic force Fjy. Such a term - of 
unclear physical origin in dense multicontact systems 
- is often introduced to ease the approach to mechan- 
ical equilibrium. Its influence will be assessed in part 
3. The viscous force is proportional to the normal rel- 
ative velocity, and the damping coefficient in the con- 
tact between grains i and j is a constant fraction ( 
(0 < C < 1) of the critical value 2( ^" t '" v ) 1 / 2 . (In 

v — => — ' \ m,i+m,j ' v 

a binary collision the normal 'restitution coefficient' 
is for ( = 1 and 1 for ( = 0). (, K N , K T , and /i 
are the same in all contacts. The motion of grains is 
calculated on solving Newton's equations. 

2.2 Numerical compression tests 

Two different types of boundary conditions are used : 
either the container walls are physical objects, with 
masses, satisfying Newton's equations (but requested 
to move in the direction perpendicular to their ori- 
entation), or periodic boundary conditions (no walls) 
are implemented. In both cases, the changes in cell 
size and shape under controlled stress involves spe- 
cific dynamical parameters which could be discussed 
in more detail. Here we will simply deem such param- 
eter choice innocuous if results are reproducible, size- 
independent and consistent. We use soil mechanics 
sign conventions for stresses and strains. Samples are 
first compressed isotropically under a constant pres- 
sure P. Once a mechanical equilibrium is reached 
under pressure P, samples are submitted to biaxial 
compression tests. The lateral stress, o\ is maintained 
equal to P, while either e 2 is increased at a con- 
stant rate e 2 (a procedure hereafter referred to as SRC, 
for strain rate controlled) or a 2 is stepwise increased 
by small fractions of P , and one waits for the next 
equilibrium configuration before changing a 2 (a SIC, 
for stress increment controlled, procedure). In the se- 
quel q denotes the ratio (a 2 — oi)/<Ji, while e 2 and 
e v = Ci + £2 — £1^2 are respectively termed 'axial' and 
'volumetric' strain, in analogy with 3D axisymmetri- 
cal triaxial tests. 

2.3 Dimensional analysis 

Rheological curves and internal sample states ob- 
tained in monotonous biaxial tests are defined in the 
macroscopic limit iV — > 00. If expressed by re- 
lations between dimensionless quantities e 2 , q, e v , 
they should depend on the friction coefficient /i and 
on ratio K T /K N , and on three other dimensionless 
parameters: k = K^/P, the stiffness parameter, 
which expresses the level of contact deformation, 7 = 

e 2 ^Jm/P, the inertia parameter, evaluating, in SRC 
(constant e 2 ) tests, the importance of dynamical ef- 
fects, and (, the damping parameter, introduced in 
paragraph ^. 1[ characterizing viscous dissipation. The 



contact coordination number is a decreasing function 
of k. The quasistatic limit is the limit of small 7. 

3 BIAXIAL COMPRESSION OF DENSE SYS- 
TEMS : RESULTS 

3.1 Preparation, initial states, procedures. 

The sample preparation procedure is well known to 
exert a strong influence on the mechanical properties 
of a granular sample as, in particular, dense or loose 
initial states respond differently (Wood 1990) to load 
increments. Moreover, experiments also showed that 
density is not sufficient to determine the behaviour 
in a triaxial test (Benahmed 2001). Numerical simu- 
lations may in principle attempt to imitate as closely 
as possible laboratory experiments. The simulations 
of such processes as deposition under gravity within 
a walled container is however difficult, as it requires 
large number of particles. Inhomogeneous states one 
obtains in such cases request samples much larger 
then a representative volume element, which is itself 
much larger than the grain size. Moreover, the transi- 
tion from an initial fluid-like configuration to a solid- 
like grain assembly is bound to be sensitive to static 
and dynamic parameters (ISilbert et al. 20011) . 

Here we focus on the slow quasistatic deforma- 
tion of certain types of granular assemblies, once 
they have been prepared in some well defined ini- 
tial state. Therefore we leave a detailed (and nec- 
essary) study of the preparation process to future 
research, and adopt a simple numerical procedure 
which provides us with homogeneous, reproducible, 
sample size -independent initial states in equilibrium 
under an isotropic pressure. The numerical procedure 
is an isotropic, monotonous compaction from an ini- 
tial gas-like configuration with a solid fraction $ of 
about 20%. To obtain a dense sample, a different, 
smaller value is attributed to the coefficient of fric- 
tion in this initial dynamic compression step. Two se- 
ries of samples are studied here. The first one - called 
series A hereafter - was prepared between solid, fric- 
tionless walls. It was observed in that case that one 
had to set p, to zero in the preparation stage if we 
were to obtain a homogeneous stress field. Simula- 
tions of series A were therefore performed starting 
from the very dense states which result from a com- 
pression without intergranular friction (Com b~e~2001l) . 
The results below, some of which were presented 
in (|Roux and Combe 20021) . were obtained with p = 
0.25 during biaxial compressions, and a rigidity level 
k = 10 5 . K T /K N was set to 1/2. Biaxial tests were 
SIC, with small q steps 5q = 10~ 3 . Each succes- 
sive mechanical equilibrium is deemed attained when 
the total force (or torque) on each grain is less than 
10~ 4 aP (resp. 10~ 4 a 2 P) and when the relative dif- 
ference between the internal overall stresses (deduced 
from non-viscous intergranular forces) and their pre- 
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Table 1 : Initial state data for series B simulations. 
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scribed values is less than 10" 4 . ( was set to high val- 
ues (near 1) and N ranged from 1024 to 4900. In the 
initial isotropic state, the solid fraction (extrapolated 
to N -> oo) is $ = 0.844 ± 0.001, all but 5.5% of 
the disks carry forces and the coordination number, 
ignoring those inactive grains, is z ~ 4.01, very close 
to the isostatic limit (|Roux 2000|) of 4 reached with 
rigid, frictionless disks in equilibrium. 

For the second series of simulations, series B, we 
used periodic boundary conditions. Samples are thus 
devoid of edge effects. They shrink homogeneously 
in the isotropic compression stage. Series B samples 
were compressed with // = 0.15, and subsequent bi- 
axial tests performed with /i = 0.5. Different stiff- 
ness levels, (k = 10 3 , 10 4 and 10 5 ) were used, with 
Kt/Kn fixed to 1, as well as different inertia param- 
eters 7 (10~ 3 , 10~ 4 , sometimes 10~ 5 ). SRC tests were 
compared to SIC ones (with 5q = 10 -2 and £ ~ 1). 
Samples of 1400 and 5600 disks were simulated. The 
initial solid fraction, due to the finite [i value during 
compression, is lower than for A samples, as well 
as the coordination number z among force-carrying 
disks. Values of $, z, and the fraction of inactive disks 
xq, for the investigated k values are given in table [TJ 
The typical aspect of q versus e 2 curves is illustrated 
on fig. [TJ for series B samples with k = 10 4 and 
7 = 10~ 4 . They are characteristic of very dense sam- 
ples, as in (IKuhn 19991) . 

3.2 Stress-strain curves and macroscopic limit. 
The increase of q with e 2 is initially quite fast, q reach- 
ing about 0.8 for e 2 < 10~ 3 . Then the deviator stress 
keeps increasing and reaches an apparent plateau for 
62 ~ 0.01. Those dense samples are markedly dilatant 
(fig. |4] below), after a very small initial contraction 
their volume steadily increases, even after q appears 
to have levelled off. The important stress fluctuations 
in those SRC tests is striking on fig. [TJ but are con- 
siderably reduced, as well as sample-to-sample dif- 
ferences, as N increases from 1400 to 5600. Dilatancy 
curves (see below) are smoother. Smooth stress-strain 
curves can thus be expected in the macroscopic limit 
N — > oo. This was more carefully checked for sim- 
ulation series A, on studying three sample sizes : on 
fig. [2] the shaded zones extend to one standard devi- 
ation on each side of the average curves, for q plot- 
ted as a function of e 2 for N = 1024 (26 samples), 
N = 3025 (10 samples), and N = 4900 (7 samples). 
Fig.[2]does indicate a systematic decrease of the fluc- 
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DEVIATORIC STRESS RATIO q 
2 SAMPLES N=1400 (THIN LINES) 
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Figure 1 : q versus axial strain e 2 in B samples of 2 differ- 
ent sizes. Fluctuations are larger for the smaller samples. 



tuation level (see inset), compatible with a regression 
as A^ -1 / 2 , just like for an average over a number of 
independent contributions (subsystems of representa- 
tive size) proportional to N. Series A samples respond 
in a similar way to deviator stresses as type B ones (al- 
though of course, due to different initial states, [i and 
k, constitutive laws will differ). The initial increase 
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Figure 2: Hashed zone (the darker the larger N) one r.m.s. 
deviation on each side of average curve for the 3 sample 
sizes indicated (series A, SIC with 5q = 10~ 3 ). Inset : its 
average width over the e 2 < 0.02 interval, versus 1/y/N, 
along with the average relative uncertainty on e v 
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of q, so fast that it cannot be distinguished from the 
axis on fig. |2l is followed by a slower variation. (Yet, 
unlike in the B case, q does not reach a maximum for 
e 2 < 0.02). 'Volumetric' strains are also qualitatively 
similar for series A and B. 

3.3 Role of parameters (, -y, k. 

The quasistatic stress-strain curve should be the same 
for SRC and SIC biaxial compressions, independent 
on £ and on 7 if it is small enough. To check this, five 
samples of series B were submitted to SRC tests with 
7 = 10~ 3 and C = 1, 7 = 10" 4 and ( = 1, 7 = 10~ 4 
and C = 0, and to SIC ones with Sq = 10~ 2 . Av- 
erage curves for q versus e 2 (fig- 13) and e v versus e 2 
(fig- S) for those 4 sets of simulations are displayed 
(and standard deviations levels indicated as on fig. [2]). 
Obviously, the value of ( does not have any apprecia- 
ble influence on the rheological curve. Intergranular 
friction is the dominating dissipation mechanism, and 
it can be checked that the differences between stresses 
evaluated with and without viscous forces differ by 
negligible amounts for all SRC tests. However, results 
are affected by the reduced rate 7, or the choice of 
an SIC procedure. A smaller 7 (according to its def- 
inition, this amounts to a slower compression, lighter 
grains or higher pressures) results in smaller devia- 
tor and dilatancy values for a given 'axial' strain. SIC 
tests, as one waits for equilibrium, are the slowest, 
and SIC curves can be regarded as an extrapolation 
of SRC ones to 7 = 0. (The occurrence of slightly 
decreasing q values in SIC tests might seem surpris- 
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Figure 3: Average q versus e 2 for conditions indicated. 
Left inset: detail of one curve with r.m.s. deviations, small 
e 2 . Right inset: averages and r.m.s. deviations for 7 = 
1(T 3 , 7 = 1(T 4 and SIC tests. 
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Figure 4: Same as fig.[3]for e v vs. e 2 , standard deviations 
shown except for uppermost (7 = 10~ 3 ) curve. 

ing, but is due to the use of real Cauchy stresses to 
draw the curve, while stresses defined in terms of ini- 
tial cell dimensions are used in the calculations). The 
effects of the stiffness parameter k are illustrated on 
fig- O It is most apparent in the initial rise of q, which 
is the faster for higher k, and the small- strain con- 
tractant regime (see inset), which develops with softer 
contacts. For smaller k, the packing appears indeed to 
be softer. The curves at larger strains display no con- 
spicuous difference between k = 10 4 and k = 10 5 , 
although the softest grains, k = 10 3 appear to with- 
stand a somewhat higher deviator stress. The dila- 
tancy - slope of — e v versus e 2 - is not affected. The 
time scale for stress fluctuation during monotonous 
tests at a given strain rate is a strongly decreasing 
function of k, hence the smoother curves on fig.[5]for 
softer contacts. The effects of the parameters on rheo- 
logical curves are related to some changes in the inter- 
nal states of the system undergoing compression. The 
effect of 7 is related to the greater distance to equilib- 
rium of systems under higher strain rate. Characteris- 
tic quantities are the average kinetic energy per parti- 
cle, e c (in units of a 2 P) and the quadratic average of 
the net force on a particle (in units of aP), f 2 . Those 
quantities tend to slowly increase with e 2 during the 
test, but typical values for e 2 = 0.01 can be cited. As 
for SIC tests, one only records equilibrium positions, 
ensuring f 2 < 10~ 5 and e c < 10~ 8 . The coordina- 
tion number z and the proportion of sliding contacts 
X s vary quickly before e 2 = 10 -3 and remain es- 
sentially constant afterwards (one has z = 3.12, on 
average, for k = 10 4 and 7 = 10~ 4 , z = 3.05 for 
k = 10 4 and 7 = 10~ 3 ). Tests with the highest 7 val- 
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ues 10 3 are, logically, the farthest from equilibrium 
(e c = 1.5 10" 5 , and f 2 = 0.01, while e c ~ 5 10~ 7 
and f 2 = 0.02 for 7 = 10~ 4 ). The change of re makes 
a significantly larger difference from 10 4 to 10 3 than 
from 10 5 to 10 4 . Unlike e c and f 2 , which essentially 
depend on 7, z and X s are sensitive to both parame- 
ters. In SIC tests (re = 10 4 ), z decreases from its initial 
value to about 3.22 (for e 2 ~ 0.01) which is consistent 
with its dependence on 7 in SRC conditions. Interme- 
diate configurations of SIC tests, remarkably, do not 
have any sliding contact: on approaching equilibrium, 
all contact forces leave the edge of the Coulomb cone. 
Upon resuming an SRC motion, very small displace- 
ments can mobilize friction and X s > is observed 
(typically X s ~ 10%, if 7 = 10" 4 and re = 10 4 , X s 
increases with 7 and with re). 

4 DIFFERENT ORIGINS OF STRAIN 

One striking aspect of the rheological curves is the 
existence of two different regimes. At small e 2 , close 
to the initial isotropic state, curves are quite smooth 
and reproducible, sample to sample fluctuations are 
very small (figs. Q] and [3]), SIC and SRC tests (what- 
ever 7 < 10 -3 ) are in perfect agreement (figs. [3] 
and HI), and re strongly affects the results (fig. [5]). Co- 
ordination numbers and friction mobilization change 
fast from initial values (table [T]) to the roughly con- 
stant ones given in paragraph 13 .31 At larger strains, 
the system is sensitive to the strain rate, much more 
than to the stiffness parameter. Fluctuations are con- 
siderably larger, and the stepwise increase of q, as 
one records the ensuing sequence of equilibria, re- 
sults in a staircase-shaped q versus e 2 curve, as on 
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Figure 5: Results for one B-sample with 3 different stiff- 
ness values, q (main plot) and e v (inset) vs. e 2 . 
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e 2 

Figure 6: Two SIC q vs. e 2 curves. Inset : initial strictly 
quasistatic regime, blown-up e scales. Results on one sam- 
ple are identical with both static and dynamic methods. 

fig- IS q increments in those SIC simulations are very 
small, 5q = 10~ 3 , so that nearly vertical segments 
on those plots correspond to many different equilib- 
rium configurations, each very close to the previous 
one. The slope of those steep parts of the curve is 
close to that of the initial, stiff rise of q, confused 
with the axis on the main plot in the figure, and visible 
in the blown-up inset. Large horizontal segments are 
due to motions between more distant configurations. 
The origin of those two different regimes is clarified 
once it is attempted to find the system response to 
small load increments by purely static means. Start- 
ing from an equilibrium configuration, it is possible 
to regard its contact structure as a given network of 
elastoplastic elements, and determine the displace- 
ments leading to the new equilibrium configuration, 
with a static method which is a discrete analog of 
elastoplastic finite element calculations in continuum 
mechanics. Such methods are seldom used (see, how- 
ever, (Kishi no et al. 200T1) ) in granular systems be- 
cause they are more complicated and less versatile 
than the usual dynamical approaches: a stiffness ma- 
trix has to be rebuilt for each different contact list, and 
calculations are limited to the range of stability of a 
given contact network. As long as the contact struc- 
ture is able to support the load, plastic strains in the 
sliding contacts remain contained by elastic strains in 
the non-sliding ones, and the static method is able to 
determine the sequence of configurations reached on, 
e.g., stepwise increasing q. This sequence is made of 
a continuous set of equilibrium states, and the system 
evolution is indeed quasistatic : we refer to such case 
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as the strictly quasistatic regime. We checked, for se- 
ries A samples, that static and dynamic calculations 
are in perfect agreement in such cases, as shown on 
fig- IB This initial regime is the stability range of the 
initial configuration. The strains are then directly due 
to contact deformation - such strains will be termed 
of type I in the sequel - and are inversely propor- 
tional to k, while results are not sensitive to 7 (the 
static method ignores completely inertia and physical 
time). This range should not be regarded as an elas- 
tic domain, as the non-linearity of the curves on fig. [6] 
(the elasticity of contacts is linear) is due to contact 
losses and also to the gradual mobilization of friction. 
On reversing the q increments, steeper slopes are ob- 
served. In the samples of fig. [6l the very steep parts 
of the staircase-shaped curves also correspond, as we 
checked, to stability intervals of some intermediate 
equilibrium configuration at higher q. Such intervals 
are separated by large strain steps, corresponding to 
rearrangements of the contact structure. Those occur 
when the accumulation of sliding contacts leads to an 
instability, and the ensuing motion is arrested by new 
contacts as interstices between neighbouring grains 
are closed. The resulting strain increments are here- 
after referred to as type II strains. Their magnitude 
is related to the width of interstices between neigh- 
bouring grains. The system evolution, in that rear- 
rangement regime, is, as shown previously, more sen- 
sitive to dynamical parameter 7. Equilibrium states 
do not form a continuum in configuration space, the 
system has to jump between two successive ones in 
a controlled deviator step test, or to flow nearby in 
a controlled strain rate test. The evolution can only 
be termed quasistatic in a wider sense if the statis- 
tical properties of trajectories in configuration space 
are independent, for slow enough motions, on dynam- 
ical parameters - which can be reasonably expected 
from the present study. The initial strictly quasi-static 
q < qi interval does not shrink, but appears rather to 
approach a finite limit (about qx = 0.8 here) as the 
sample size increases. Stress-strain curves depend on 
K T /K N within this range, but, interestingly, qx does 
not (ICombe 20011) . In the rearrangement regime, in 
order to approach a smooth curve in the macroscopic 
limit (see fig. EK it is necessary that the sizes of both 
the steep and the flat parts of the 'staircases' shrink 
to zero as the sample size increases. Type I and type 
II strains have very different amplitudes in A samples 
with k = 10 5 and iV < 4900. It might in fact be ex- 
pected that this clearcut distinction will get blurred 
at smaller stiffness parameter k (whence larger type 
I strains) or larger N (as smaller type II strain incre- 
ments can close contacts), and that the transition at 
qx will be fuzzier. Nevertheless, the system proper- 
ties do strongly differ for q < qi and q > qx, in two 
important respects. First, the slope of the stress-strain 
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Figure 7: 'Creep tests', dots on main plot showing ini- 
tial and final (equilibrium) states. Effect of resuming com- 
pression SRC way shown as thick lines. Inset: creep tests 
within strictly quasistatic range. 

curve relates directly to the elasticity of the contacts 
in the type I strain dominated, strictly quasistatic case. 
The tangent at the origin on fig. [6] (smaller plot) is the 
Young modulus of the packing. Second, the ampli- 
tude of fluctuations, the distance to mechanical equi- 
librium, and the sensitivity to perturbations are much 
stronger in the rearrangement (type II strain domi- 
nated) regime. This is further illustrated by the fol- 
lowing 'creep experiment' : in a strain-rate controlled 
biaxial compression, at some arbitrary instant, shift to 
stress-controlled conditions and keep q constant, until 
an equilibrium configuration is reached. Typical re- 
sults of such tests are shown on fig. [7J As could be 
expected, much larger strain variations are observed 
during periods of creep in the rearrangement regime, 
as the initial states are farther from equilibrium. On 
resuming the constant strain rate test, the initial part 
of the curve is very steep, which is characteristic of 
a 'strictly quasistatic' interval. From an equilibrium 
state (devoid of sliding contacts), friction has to be 
mobilized again to produce the instabilities of the re- 
arrangement regime. The dilatancy within those creep 
intervals is similar to the SRC one. 

The 'creep tests' reveal different behaviours in the 
two deformation regimes in SRC tests. One might 
also probe the sensitivity to perturbations of inter- 
mediate equilibrium states obtained in SIC tests. We 
repeatedly applied on the grains constant external 
forces, each force component being randomly chosen 
between — f and f (/ is a small fraction of aP), un- 
til new, perturbed equilibria were reached. Such ran- 
dom load increments always tend to produce strains 
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Figure 8: Effect of repeated random load (/o/aP = 
2.5 10 -3 ) applied in states shown as big dots on the stress- 
strain curve in the inset: increments of e v vs. increments of 
62 on blown-up (by 10 5 ) scale. The response of state A is 
concentrated near the origin, only the response of state B 
(q = 0.94) is visible on this scale. Dotted lines: SIC and 
SRC dilatancy curves near point B, same sample. 

in the same direction, as illustrated on fig. [HI Applied 
when q = 0.5 within the strictly quasistatic range, 
such perturbations entail very small strain increments 
(hardly visible near the origin of the plot). Applied 
when q = 0.94 as equilibrium states are much more 
unstable, they produce the series of strain increments 
plotted as connected dots, which tend to accumu- 
late proportionnally, hence the nearly straight line, the 
slope of which is comparable to the dilatancy. The re- 
peated application of small random perturbations thus 
entails some 'creep' phenomenon. 

5 COMPARISONS WITH EXPERIMENTS 

In spite of the many differences between the numeri- 
cal models and the materials studied in the laboratory, 
such as sand, or even glass beads, some features of the 
simulation results can be compared in a qualitative or 
semi-quantitative way to experimental ones. 

First, parameters n and 7 should be used to ob- 
tain robust estimations of orders of magnitude. In 3 
dimensions, k should be defined as K N /(aP) in the 
case of linear elasticity in the contacts, n measures 
the normal elastic deflection in a contact, relatively to 
the grain diameter a, due to the typical contact force 
Pa 2 . In a Hertzian contact between spheres of diam- 
eter a, it is easy to show that k should be defined as 
(E/P) 2 / 3 , where E is the Young modulus of the grain 
material. This gives k ~ 6000 for glass beads un- 



der P = 10 5 Pa. (In 3D simulations, we could check 
that, given these definitions, the effect of k on the 
coordination number was similar to the 2D case, see 
also (IMakse et al. 20001) ). 'Real' materials with Hertz 
contacts under P = 10 5 Pa are rather on the rigid side, 
but not quite in the rigid limit. Other contact laws 
might lead to even smaller stiffness parameters (e.g., 
k ~ (E/P) 1 / 2 if F N oc Eh 2 , as for cone-shaped as- 
perities). 

An appropriate 3D definition of 7 is ey^§? (\f^ 
is the time for a grain accelerated from rest by the 
typical force a 2 P to move on distance a/2). Substi- 
tuting typical values - a fraction of millimetre for a, 
10~ 5 s _1 for e - this yields 7 values as small as 10~ 9 
or 10 -10 . As calculations over e = 2% strain intervals 
with 7 = 10~ 5 still require several days of c.p.u. time 
with 5000 stiff grains, real time scales of quasistatic 
laboratory tests are still beyond the reach of discrete 
numerical simulations. 7 dependences of numerical 
results can however be extrapolated to smaller values. 

Although it is tempting, in view of the 
results illustrated on fi g. [7] to refer to 

creep experiments (|Matsushita et al. 1 999; 

IDi Benedetto and Tatsuoka 19971) , as the aspects 
of the stress-strain curves are quite similar in several 
respects, this difference of time scales precludes 
a direct comparison. Moreover, the experimen- 
tal q-e curves do not depend on strain rate if it 
is constant (this corresponds to much smaller 7 
values than simulations), and the creep defor- 
mation is extremely slow, often logarithmic in 



time (Di Prisco and Imposimato 1997). Unlike in the 
numerical case, it does not appear to stop as some 
equilibrium is reached. It might well be relevant, 
however, to discuss such experiments in terms of the 
sensitivity of the system to perturbations, which is 
likely to depend on whether contact networks resist 
load increments (strictly quasistatic case) or are prone 
to instabilities (rearrangement regime). The numeri- 
cal tests discussed in connection with fig. [8] suggest a 
possible microscopic origin of such slow evolutions 
over long times: a small noise level, always present in 
experiments, could entail an accumulation of strain. 
Aging and creep phenomena can also be physically 
expected within one contact. Numerical simulations 
(devoid of such features) might help assessing the 
collective aspects of the packing response. 

Our simulations can also be likened to ex- 
perimental observations about the very small 
strain elastic behaviour of granular sys- 
tems (|Di Be nedett o et al. 19 99). Recent develop- 
ments of precision apparati enabled measurements of 
strains in the 10 -5 range. To obtain elastic moduli, 
small stress cycles are superimposed on a constant 
loading, producing cyclic strains on top of a sys- 
tematic drift which, on increasing the number of 
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cycles, gradually slows down and becomes analogous 
to the one observed in creep tests. The average 
slope of a cycle on a stress- strain plot, once the 
effect of the drift is negligible, can be interpreted 
as an elastic modulus (there remaining some small 
dissipation). Those small strain increment elastic 
constants agree with the ones deduced from acoustic 
wave velocities. From our simulations, it transpires 
that the incremental stress-strain dependence might 
express a genuinely elastic behaviour (supplemented 
by some plastic dissipation which vanishes in the 
limit of small stress increments) in the strictly 
quasistatic regime. Elastic moduli are then related 
to the stiffness of the contacts. The width of strictly 
quasistatic strain intervals are of the order of qi/n - 
qi being their width in terms of stress ratio. Taking 
into account that q± is exceptionally large for the 
initial small-strain regime if our extremely dense and 
well coordinated systems, and the value k ~ 6000 
estimated above for glass beads, one does obtain the 
right order of magnitude (< 10~ 4 ) for the very small 
strain elastic domain. Moreover, the procedure by 
which these moduli are measured can be interpreted 
as the preparation, either left to random perturbations 
or forced by cyclic load increments, of a better 
stabilized state for which the contact network is able 
to resist small, but finite stress increments (just like 
the stiffly responding equilibrium states of fig. [7]). 

6 CONCLUSIONS AND PERSPECTIVES 

Despite their limitations (due to the simplicity of the 
contact model, and the inaccessibility of long time 
scales), the numerical simulation results presented 
here enable some investigation of the microscopic ori- 
gins of many features of experimentally observed be- 
haviours. The definition of reduced dimensionles pa- 
rameters (k and 7) provides a framework in which 
many experimental and numerical studies can be dis- 
cussed in common terms. Due to the small size of 
numerical samples, constitutive laws have to be ap- 
proached via statistical analyses. Most importantly, 
the distinction between two different origins of strain 
and two deformation regimes allows us some inter- 
pretations of very small strain (tangential) elasticity 
and slow deformation (creep) under constant load, in 
terms of the system sensitivity to perturbations. 

This work should be pursued in three directions. 
First, it is desirable to extend the existing approach 
to more 'realistic' models, so that more quantitative 
comparisons with experiments will be possible (our 
3D results on spheres - an obvious step in this di- 
rection, were not presented here for lack of space). 
Secondly, the importance of the initial state and of 
the sample preparation procedure calls for systematic 
studies (unlike for quasistatic monotonous compres- 
sion tests, experimental knowledge is not expressed 



as well established laws for such processes). And, fi- 
nally, the joint use of dynamic and static methods, 
which agree remarkably in strictly quasistatic do- 
mains (fig. [6]) opens avenues to explore fundamental 
issues, such as elastoplastic contact network stability 
and rearrangements, in some microscopic detail. 
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